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Abstract 

Designs conditions for marine structures are typically informed by threshold-based 
extreme value analyses of oceanographic variables, in which excesses of a high threshold 
are modelled by a generalized Pareto (GP) distribution. Too low a threshold leads 
to bias from model mis-specification; raising the threshold increases the variance of 
estimators: a bias-variance trade-off. Many existing threshold selection methods do not 
address this trade-off directly, but rather aim to select the lowest threshold above which 
the GP model is judged to hold approximately. In this paper Bayesian cross-validation 
is used to address the trade-off by comparing thresholds based on predictive ability at 
extreme levels. Extremal inferences can be sensitive to the choice of a single threshold. 

We use Bayesian model-averaging to combine inferences from many thresholds, thereby 
reducing sensitivity to the choice of a single threshold. The methodology is applied to 
significant wave height datasets from the northern North Sea and the Gulf of Mexico. 

Keywords. Cross-validation; Extreme value theory; Generalized Pareto distribution; Pre¬ 
dictive inference; Threshold 

1 Introduction 

Ocean and coastal structures, including breakwaters, ships and oil and gas producing facilities 
are designed to withstand extreme environmental conditions. Marine engineering design codes 
stipulate that estimated failure probabilities of offshore structures, associated with one or more 
return periods, should not exceed specified values. To characterize the environmental loading 
on an offshore structure, return values for winds, waves and ocean currents corresponding to 
a return period of typically 100 years, but sometimes to 1,000 and 10,000 years are required. 
The severity of waves in a storm is quantified using significant wave height. Extreme value 
analyses of measured and hindcast samples of significant wave height are undertaken to derive 
environmental design criteria, typically by fitting a GP distribution to excesses of a high 
threshold. The selection of appropriate threshold(s) is important because inferences can be 
sensitive to threshold. 
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1.1 Storm peak significant wave height datasets 

The focus of this paper is the analysis of two sequences of hindcasts of storm peak significant 
wave height, shown in Figure 1. Significant wave height (H s ) is a measure of sea surface 
roughness. It is defined as four times the standard deviation of the surface elevation of a 
sea state , the ocean surface observed for a certain period of time (3 hours for our datasets). 
Cardone et al. (2014) gives the largest H s value ever generated by a hindcast model as 18.33m 
and the largest value ever measured in the ocean as 20.63m. Hindcasts are samples from 
physical models of the ocean environment, calibrated to observations of pressure, wind and 
wave fields. 

For each of the datasets raw data have been declustered, using a procedure described in 
Ewans and Jonathan (2008), to isolate cluster maxima (storm peaks) that can reasonably 
be treated as being mutually independent. We also assume that storm peaks are sampled 
from a common distribution. Even in this simplest of situations practitioners have difficulty in 
selecting appropriate thresholds. The first dataset, from an unnamed location in the northern 
North Sea, contains 628 storm peaks from October 1964 to March 1995, but restricted to the 
period October to March within each year. The other dataset contains 315 storm peaks from 
September 1900 to September 2005. 
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Figure 1: Storm peak significant wave height hindcast datasets. Left: North Sea data (628 
observations). Right: Gulf of Mexico data (315 observations). Top: time series plots. Bottom: 
histograms. The upper axis scales give the sample quantile levels. 


In the northern North Sea the main fetches are the Norwegian Sea to the North, the 
Atlantic Ocean to the west, and the North Sea to the south. Extreme sea states from the 
directions of Scandinavia to the east and the British Isles to the southwest are not possible, 
owing to the shielding effects of these land masses. At the location under consideration, the 
most extreme sea states are associated with storms from the Atlantic Ocean (Ewans and 
Jonathan, 2008). With up to several tens of storms impacting the North Sea each winter, the 
number of events for analysis is typically larger than for locations in regions such as the Gulf 
of Mexico, where hurricanes produce the most severe sea states. Most hurricanes originate in 
the Atlantic Ocean between June and November and propagate west to northwest into the 
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Gulf producing the largest sea states with dominant directions from the southeast to east 
directions. It is expected that there is greater potential for very stormy sea conditions in the 
Gulf of Mexico than in the northern North Sea and therefore that the extremal behaviour is 
different in these two locations. 


1.2 Extreme value threshold selection 


Extreme value theory provides asymptotic justification for a particular family of models for 
excesses of a high threshold. Let X 1; X 2 , ■ ■ ■ X n be a sequence of independent and identically 
distributed random variables, with common distribution function H, and u n be a threshold, 
increasing with n. Pickands (1975) showed that if there is a non-degenerate limiting distribu¬ 
tion for appropriately linearly rescaled excesses of u n then this limit is a Generalized Pareto 
(GP) distribution. In practice, a suitably high threshold u is chosen empirically. Given that 
there is an exceedance of u, the excess Y — X — u is modelled by a GP(oy,£) distribution, 
with positive threshold-dependent scale parameter <r u , shape parameter £ and distribution 
function 

i - (i + £y/^)+ i/? , i + o, 

1 - exp(—y/ (T u ), £ = 0, 



where y > 0, x + — max(x, 0). The £ = 0 case is defined in the limit as £ —» 0. When £ < 0 
the distribution of X has a finite upper endpoint of u — <r u /£; otherwise, X is unbounded 
above. The frequency with which the threshold is exceeded also matters. Linder the current 
assumptions the number of exceedances u has a binomial(n, p u ) distribution, where p u = 
P(X > u), giving a BGP(p n , a u , £) model (Coles, 2001, chapter 4) 

Many threshold diagnostic procedures have been proposed: Scarrott and MacDonald (2012) 
provides a review. Broad categories of methods include: assessing stability of model parameter 
estimates with threshold (Drees et ah, 2000; Wadsworth, 2015); goodness-of-ht tests (Davison 
and Smith, 1990; Dupuis, 1998); approaches that minimize the asymptotic mean-squared error 
of estimators of £ or of extreme quantiles, under particular assumptions about the form of the 
upper tail of H (Hall and Welsh, 1985; Hall, 1990; Ferreira et ah, 2003; Beirlant et ah, 2004); 
specifying a model for (some or all) data below the threshold (Wong and Li, 2010; MacDonald 
et ah, 2011; Wadsworth and Tawn, 2012). In the latter category, the threshold above which 
the GP model is assumed to hold is treated as a model parameter and threshold uncertainty 
is incorporated by averaging inferences over a posterior distribution of model parameters. In 
contrast, in a single threshold approach threshold level is viewed as a tuning parameter, whose 
value is selected prior to the main analysis and is treated as fixed and known when subsequent 
inferences are made. 

Single threshold selection involves a bias-variance trade-off (Scarrott and MacDonald, 
2012): the lower the threshold the greater the estimation bias clue to model misspecihcation; 
the higher the threshold the greater the estimation uncertainty. Many existing approaches 
do not address the trade-off directly, but rather examine sensitivity of inferences to thresh¬ 
old and/or aim to select the lowest threshold above which the GP model appears to hold 
approximately. We seek to deal with the bias-variance trade-off based on the main purpose 
of the modelling, i.e. prediction of future extremal behaviour. We make use of a data-driven 
method commonly used for such purposes: cross-validation (Stone, 1974). We consider only 
the simplest of modelling situations, i.e. where observations are treated as independent and 
identicaly distributed. However, selecting the threshold level is a fundamental issue for all 
threshold-based extreme value analyses and we anticipate that our general approach can have 
much wider applicability. 

We illustrate some of the issues involved in threshold selection by applying to the significant 
wave height datasets two approaches that assess parameter stability. In the top row of Figure 
2 maximum likelihood (ML) estimates £ of £ are plotted against threshold. The aim is to 
choose the lowest threshold above which £ is approximately constant in threshold, taking into 
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account sampling variability summarized by the confidence intervals. It is not possible to 
make a definitive choice and different viewers may choose rather different thresholds. In both 
of these plots our eyes are drawn to the approximate stability of the estimates at around the 
70% sample quantile. One could argue for lower thresholds, to incur some bias in return for 
reduced variance, but it is not possible to assess this objectively from these plots. In practice, 
it is common not to consider thresholds below the apparent mode of the data because the GP 
distribution has it’s mode at the origin. For example, based on the histogram of the Gulf of 
Mexico data in Figure 1 one might judge a threshold below the 25% sample quantile to be 
unrealistic. However, we will consider such thresholds because it is interesting to see to what 
extent the bias expected is offset by a gain in precision. 

The inherent subjectivity of this approach, and other issues such as the strong dependence 
between estimates of £ based on different thresholds, motivated more formal assessments of 
parameter stability (Wadsworth and Tawn, 2012; Northrop and Coleman, 2014; Wadsworth, 
2015). The plots in the bottom row of Figure 2 are based on Northrop and Coleman (2014). 
A subasymptotic piecewise constant model (Wadsworth and Tawn, 2012) is used in which the 
value of £ may change at each of a set of thresholds, here set at the 0%, 5%, ..., 95% sample 
quantiles. For a given threshold the null hypothesis that the shape parameter is constant 
from this threshold upwards is tested. In the plots p-values from this test are plotted against 
threshold. Although these plots address many of the inadequacies of the parameter estimate 
stability plots subjectivity remains because one must decide how to make use of the p-values. 
One could prespecify a size, e.g. 0.05, for the tests or take a more informal approach by 
looking for a sharp increase in p-value. For the North Sea data the former would suggest a 
very low threshold and the latter a threshold in the region of the 70% sample quantile. For 
the Gulf of Mexico data respective thresholds close to the 10% and 55% sample quantiles are 
indicated. 
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Figure 2: Threshold diagnostic plots for the storm peak significant wave height hindcast 
datasets. Left: North Sea data. Right: Gulf of Mexico data. Top: parameter stability plots 
for MLEs of £, with 95% pointwise profile likelihood-based confidence intervals. Bottom: 
p-values associated with a test of constant shape parameter against the lowest threshold 
considered. The upper axis scales give the level of the threshold in metres. 
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An argument against selecting a single threshold is that this ignores uncertainty concerning 
the choice of this threshold. As mentioned above, one way to account for this uncertainty 
is to embed a threshold parameter within a model. We use an approach based on Bayesian 
Model Averaging (BMA), on which Hoeting et al. (1999) provide a review. Sabourin et al. 
(2013) have recently used a similar approach to combine inferences from different multivariate 
extreme value models. We treat different thresholds as providing competing models for the 
data. Predictions of extremal behaviour are averaged over these models, with individual 
models weighted in proportion to the extent to which they are supported by the data. There 
is empirical and theoretical evidence (Hoeting et ah, 1999) that averaging inferences in this 
way results in better average predictive ability than provided by any single model. 

For the most part we work in a Bayesian framework because prediction is handled natu¬ 
rally and regularity conditions required for making inferences using ML (Smith, 1985) and 
probability-weighted moments (PWM) (Hosking and Wallis, 1987), namely £ > —1/2 and 
£ < 1/2 respectively, can be relaxed. This requires a prior distribution to be specified for the 
parameters of the BGP model. Initially, we consider three different ‘reference’ prior distri¬ 
butions, in the general sense of priors constructed using formal rules (Kass and Wasserman, 
1996). Such priors can be useful when information provided by the data is much stronger 
than prior information from other sources. This is more likely to be the case for a low thresh¬ 
old than for a high threshold. We use simulation to assess the utility of these priors for our 
purpose, that is, making predictive statements about future extreme observations and use 
the results to formulate an improved prior. For high thresholds, when the data are likely to 
provide only weak information, it may be important to incorporate at least some basic prior 
information in order to avoid making physically unrealistic inferences. 

In section 2 we use Bayesian cross-validation to estimate a measure of threshold perfor¬ 
mance for use in the selection of a single threshold. In section 3 we consider how to use this 
measure to combine inferences over many thresholds. Sections 2.1 and 2.2 describe the cross- 
validation procedure and its role in selecting a single threshold. In section 2.3 we discuss two 
related formulations of the objective of an extreme value analysis and, in section 2.4, we use 
one of these formulations in a simulation study to inform the choice of prior distribution for 
GP parameters. Another simulation study, in section 3.1, compares the strategies of choosing 
a single ‘best’ threshold and averaging of many thresholds. In sections 2.5 and 3.2 we use our 
methodology to make inferences about extreme significant wave heights in the North Sea and 
in the Gulf of Mexico. In section 3.3 we illustrate how one could incorporate prior information 
about the GP shape parameter £ to avoid physically unrealistic inferences. 

2 Single threshold selection 

We use a Bayesian implementation of leave-one-out cross-validation to compare the predictive 
ability of BGP inferences based on different thresholds. We take a predictive approach, aver¬ 
aging inferences over the posterior distribution of parameters, to reflect differing parameter 
uncertainties across thresholds: uncertainty in GP parameters will tend to increase as the 
threshold is raised. A point estimate of GP model parameters can give a zero likelihood for a 
validation observation: this occurs if £ < 0 and this observation is greater than the estimated 
upper endpoint u — a u /£- In this event an estimative approach would effectively rule out the 
threshold u. Accounting for parameter uncertainty alleviates this problem by giving weight 
to parameter values other than a particular point estimate. 

A naive implementation of leave-one-out cross-validation is computationally intensive. To 
avoid excessive computation we use importance sampling to estimate cross-validation predic¬ 
tive densities based on Bayesian inferences from the entire dataset. One could use a similar 
strategy in a frequentist approximation to predictive inference based on large sample the¬ 
ory or bootstrapping (Young and Smith, 2005, chapter 10). However, large sample results 
may provide poor approximations for high thresholds (small number of excesses) and the GP 
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observed information is known to have poor finite-sample properties (Siiveges and Davison, 
2010). Bootstrapping, of ML or PWM estimates, increases computation time further and is 
subject to the regularity conditions mentioned in the introduction. 


2.1 Assessing threshold performance using cross-validation 

Suppose that x = (xi,...,x n ) is a random sample of raw (unthresholded) data from H. 
Without loss of generality we assume that x\ < • • • < x n . Consider a training threshold 
u. A BGP(p n , cr u , £) model is used at threshold u, where p u = P(X > u ) and (cr u ,£) are 
the parameters of the GP model for excesses of u. Let 9 = (p u , v ui £) and tt(9) be a prior 
density for 6. Let x s denote a subset of x , possibly equal to x. The posterior density 
7 t u (9 | x s ) oc L(9;x s ,u)ir(9), where 

L(9; x s , u) = JJ f u (xi\9), 

i:Xi£x s 

fu{xi | 9) = (1 - Pu) I{Xi ^ u) {p u g{xi ~ u: a u ,0} I{Xi>u) , 

I(x) = 1 if x is true and I(x) = 0 otherwise, and 

, -1A + ^v (1+1/€) 

g{x](r u ,^) = a u 1 + — 

V a uJ + 

is the density of a GP(<r n ,£) distribution. 

We quantify the ability of BGP inferences based on threshold u to predict (out-of-sample) at 
extreme levels. For this purpose we introduce a validation threshold v ^ u. If 1 +£(v—u)/a u > 
0 then a BGP(p u , a u , £) model at threshold u implies a BGP(p„, a v , £) model at threshold v, 
where a v = a u + £(v — u) and p v = P(X > v) = (1 + £{v — u)/a u )~ 1 ^p u . Otherwise, p v — 0 
and excesses of v are impossible. For a particular value of v we wish to compare the predictive 
ability of the implied BGP(p t ,, a v , £) model across a range of values of u. 

We employ a leave-one-out cross-validation scheme in which aj( r ) = {xi,i ^ r} forms the 
training data and x r the validation data. The cross-validation predictive densities at validation 
threshold v, based on a training threshold u, are given by 


f v (x r I X( r ),u) = I f v (x r \ 9, £c (r) ) 7t u (9 \ *( r )) d 9, r = l,...,n. (2) 

Suppose that the {xj} are conditionally independent given 9. If p v > 0 then 

fv(x r | 9, X(r)) = f v (x r | 9) = (1 - p v ) I{ ~ Xr ^ v) {p v g(x r - v; a v , ^)} / (^>D _ (3) 

If p v — 0 then f v (x r \ 9 , x^) — I(x r ^ v). Suppose that we have a sample 0'p , j = 1 ,pm 
from the posterior 7 t u {0 \ x^). Then a Monte Carlo estimator of f v (x r \ x^,u) based on (2) 
is given by 



m 


5 ~2fv(x r \ X(r))- 


3 = 1 


(4) 


Evaluation of estimator (4), for r = 1,... , n, is computationally intensive because it involves 
generating samples from n different posterior distributions. To reduce computation time we 
use an importance sampling estimator (Gelfand, 1996; Gelfand and Dey, 1994) that enables 
estimation of f v (x r \ x^,u), for r = 1,..., n — 1, using a single sample only. We rewrite (2) 
as 


fv(.Xr | X( r ), tlj 


I f v (x r | 9, ® (r) ) q r {9 ) h(9) d 9, r = 1,..., n, 


(5) 
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where q r (6) = tt u {6 j X( r ))/h(6) and h{9) is a density whose support must include that of 
7 t u {6 | £C( r )). In the current context a common choice is tt u (6 \ x ) (Gelfand and Dey, 1994, 
page 511). However, the support of n u (9 \ x): £ > —<J u /(x n — u), does not contain that of 
n u (9 | £C( n )): £ > — cr u /(x n -i — u), since x n > x n -\. Therefore we use h{6) = n u {9 \ x) only 
for r % n. 

Suppose that we have a sample 9j,j = 1 , m from the posterior n u (9 \ x). For r = 
1 ,... 7 n — 1 we use the importance sampling ratio estimator 



EJLl fv( x r | 0j)?r(0j) 

EJLl Qr(0j) ’ 

EJlii/Ou- I Oj)/fu(x r 10j) 
' E7=l VWr I 


( 6 ) 

(7) 


where q r {9) = n u {9 \ x^)/i t u (9 \ x) oc 1/f u {x r \ 9). If we also have a sample 9^\j = 1 ,m 

from the posterior n u (9 j *(„)) then f v (x n \ x^,u) = (1/m) Y^=i fv{?n \ 9 j n) ). 

We use 

n 

T v {u ) = y,log f v (x r I *( r ),M) (8) 

r=l 

as a measure of predictive performance at validation threshold v when using training threshold 
u. 


2.2 Comparing training thresholds 

Consider k training thresholds U\ < ■ ■ ■ < Uk , resulting in estimates T v (u\),. . . , T v (uk). Up to 
an additive constant, T v (u) provides an estimate of the negated Kullback-Leibler divergence 
between the BGP model at validation threshold v and the true density (see, for example, Sil¬ 
verman (1986, page 53)). Thus, u* = argma x u T v (u) has the property that, of the thresholds 
considered, it has the smallest estimated Kullback-Leibler divergence. 

Some inputs are required: u\,... ,Uk, v and 7r(9). Choosing a set ui,... ,Uk of plausible 
thresholds that span the range over which the bias-variance trade-off is occurring is the 
starting point for many threshold selection methods. An initial graphical diagnostic, such as 
a parameter stability plot, can assist this choice. In particular, Uk should not be so high that 
little information is provided about GP parameters. There is no definitive rule for limiting 
Uk but Jonathan and Ewans (2013) suggest that there should be no fewer than 50 threshold 
excesses. Applying this rule would restrict Uk to be no higher than the 84% and 92% sample 
quantiles for the Gulf of Mexico and North Sea datasets respectively, but later we will use Uk 
that break the rule and examine the consequences. 

We need v ^ Uk, but the larger v is the fewer excesses of v there are and the smaller 
the information from data thresholded at v. Consider two validation thresholds: V\ = Uk 
and V r 2 > Uk . If we use ip we lose validation information: if v\ < x r ^ ip then in (3) x r is 
censored rather than entering into the GP part of the predictive density; and gain nothing: 
the prediction of x r > ip is unaffected by the choice of iq or ip because p vi g(x — iq; a vi ; £) = 
p v 2 g(x — ip; cr,., 2 ; £). Therefore, we should use v = Uk- 

In section 2.4 we compare predictive properties of three ‘reference’ priors for GP param¬ 
eters. Such priors are intended for use when substantial prior information is not available 
and it is anticipated that information provided by the data will dominate the posterior dis¬ 
tribution (O’Hagan, 2006). For high thresholds this may not be the case, with the possible 
consequence that a diffuse posterior places non-negligible posterior probability on unrealisti¬ 
cally large values of £ and produces unrealistic extrapolation to long future time horizons. In 
this event one needs to incorporate more information, perhaps by using a more considered 
prior distribution, and/or limit the time horizon of interest. We return to this issue in section 
3.3. 
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2.3 Prediction of extreme observations 


In an extreme value analysis the main focus is often the estimation of extreme quantiles 
called return levels. Let denote the largest value observed over a time horizon of N years. 
The N- year return level z(N) is defined as the value exceeded by an annual maximum Mi 
with probability 1/N. In off-shore engineering design criteria are usually expressed in terms 
of return levels, for values of N such as 100, 1000, 10,000. A related approach defines the 
quantity of interest as the random variable Mn, rather than particular quantiles of Mi. Linder 
a BGP(p u , cr u , £) model, for z > u, 


F(z\ D) = P(X < z) = 1 - p„ jl + « (^) } ^ ■ 

Then z(N) = z(N;9) satisfies F(z(N); 9) ny = 1 — 1/N, where n y is the mean number of 
observations per year. Similarly, for z > u, P(M ^ z) — F(z; 9) ,lyN . For large N (N = 100 
is sufficient), z(N) is approximately equal to the 37% quantile of the distribution of Mn 
(Cox et ah, 2002). In an estimative approach, based on a point estimate of 6, the value of 
z(N) is below the median of Mn- A common interpretation of z(N) is the level exceeded 
on average once every N years. However, for large N (again N = 100 is sufficient) and 
under an assumption of independence at extreme levels, z(N) is exceeded 0, 1, 2, 3, 4 times 
with respective approximate probabilities of 37%, 37%, 18%, 6% and 1.5%. It may be more 
instructive to examine directly the distribution of Mn, rather than very extreme quantiles of 
the annual maximum M\. 

The relationship between these two approaches is less clear under a predictive approach, 
in which posterior uncertainty about 6 in incorporated into the calculations. The N-year 
(posterior) predictive return level zp(N) is the solution of 


P(Mi ^ z P (N) \x)= F(z P (N)]0) ny vr (9 \ x) 6.9 = 1 - 1/N, 


and the predictive distribution function of M N is given by 


P(M n ^z\x)= I F(z;9) nyN n{9 \ x) d 9 


(9) 


As noted by Smith (2003, section 1.3), accounting for parameter uncertainty tends to lead 
to larger estimated probabilities of extreme events, that is, zp(N) tends to be greater than 
an estimate 'z(N) based on, for example, the MLE 9. The strong non-linearity of F(z; 9) Uy 
for large z, and the fact that it is bound above by 1, mean that averages of F(z ; 9) ny over 
areas of the parameter space relating to the extreme upper tail of M\ tend to be smaller than 
point values near the centre of such areas. This phenomenon is less critical when working 
with the distribution of Mn because now central quantiles of Mn also have relevance, not just 
particular extreme tail probabilities. Numerical results in section 2.5 (Figure 7) show that 
Zp(N) can be rather greater than the median of the predictive distribution of Mjv, particularly 
when posterior uncertainty about 9 is large. 

For a given value of N , we estimate P(M N ^ z \ x) using the sample 9 3 , j = 1,... pm from 
the posterior density n(9 \ x) to give 


^ 1 
P(M n ^z I *) = — V F(z- 9A nyN . 
m 

3 =i 

The solution zp{N ) of P(Mi % z'p(N) \ x) — 1 — 1/N provides an estimate of zp{N). 


( 10 ) 



2.4 Simulation study 1: priors for GP parameters 

We compare approaches for predicting future extreme observations: a predictive approach 
using different prior distributions and an estimative approach using the MLE. We use Jeffreys’ 
prior p u ~ beta(l/2,1/2) for p u , so that p u \ x ~ beta(n u + 1/2, n — n u + 1/2), where n u 
is the number of threshold excesses. Initially we consider three prior distributions for GP 
parameters: a Jeffreys’ prior 

+ 0“ 1 ( 1 + 2 0 ~ 1/2 , > 0 , f > - 1 / 2 ; ( 11 ) 

a maximal data information (MDI) prior 

km((7u,0 oc a~ 1 e~ ( ^ +1) a u > 0 , f >- 1 , ( 12 ) 

truncated from ^ 6 R to ( ^ — 1; and a flat prior 

n F ((Tu,0 oc a” 1 , a u > 0, f G M. (13) 

Motivated by Endings presented later in this section we generalize (12) to an MDI(a) prior: 

ita^, G; a ) cr“ 1 ae _a ^ +1 - ) cr u > 0, £ ^ — 1, a > 0. (14) 

These priors are improper. Let n u be the number of threshold excesses. Castellanosa and 
Cabras (2007) show that the Jeffreys’ prior yields a proper posterior for n u ^ 1 and Northrop 
and Attalides (2014) show that under the flat prior a sufficient condition for posterior propriety 
is n u ^ 3. Northrop and Attalides (2014) also show that for any sample size, if, and only if, 
£ is bounded below a priori, the MDI prior, and the generalized MDI prior, yield a proper 
posterior. At the particular bound of —1 used in (12) the GP distribution reduces to a uniform 
distribution on (0, a) and corresponds to a change in the behaviour of the GP density: for 
£ < — 1, this density increases without limit as it approaches its mode at the upper end point 
—<t m /£, behaviour not expected in extreme value analyses. Figure 3 compares the Jeffreys’, 
MDI and generalized MDI prior (for a = 0.6) as functions of £. The Jeffreys’ prior (11) is 


>% 
c n 
c: 
<D 
T3 

O 

Q. 


— Jeffreys' 
-- MDI (1) 
- MDI (0.6) 
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Figure 3: Jeffreys’, truncated MDI and generalized MDI priors as functions of £. 

unbounded as £ — 1/2. If there are small numbers of threshold excesses this can result in a 

bimodal posterior distribution, with one mode at £ = —1/2. In this simulation study we also 
find that the Jeffreys’ prior results in poorer predictive performance than the truncated MDI 
and flat priors. 

Let Z new be a future N -year maximum, sampled from a distribution with distribution 
function F(z; 6) nyN . If the predictive distribution function (9) is the same as that of Z new then 
P(Mn ^ Z new | x) has a U(0,1) distribution. In practice this can only hold approximately: 
the closeness of the approximation under repeated sampling provides a basis for comparing 
different prior distributions. Performance of an estimative approach based on the MLE 6 can 
be assessed using F(Z new ; 9) nyN . For a given prior distribution and given values of N,n y and 
n, the simulation scheme is: 
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1. simulate a dataset cc S j m of n independent observations from a BGP(p n , <r u , £) model and 
then a sample 0 3 , j = 1,pm from the posterior tt(9 | a3 sim ); 

2. simulate an observation 2 new from the distribution of Mjv, that is, max(Xi,... , AWJ, 
where N u ~ bin (n y N,p u ) and X t GP(cr u , £), i — 1,..., N u ; 

3. use (10) to evaluate P(Mn ^ ^ new | x). 

Steps 1. to 3. are repeated 10,000 times, providing a putative sample of size 10,000 from a 
U(0,1) distribution. In the estimative approach step 3 is replaced by evaluation of F(z new ] 6) n y N 
Here, and throughout this paper, we produce samples of size m from the posterior distribution 
7 t{6 | x) using the generalized ratio-of-uniforms method of Wakefield et ah (1991), following 
their suggested strategy of relocating the mode of n(6 \ x ) to the origin and setting a tuning 
parameter r to 1/2. In the simulation studies we use m = 1,000 and when analyzing real 
data we use m = 10, 000. 

We assess the closeness of the U(0,1) approximation graphically (Geweke and Amisano, 
2010), comparing the proportion of simulated values in each U(0,1) decile to the null value 
of 0.1. To aid the assessment of departures from this value we superimpose approximate 
pointwise 95% tolerance intervals based on number of points within each decile having a 
bin(10,000, 0.1) distribution, i.e. 0.1 ± 1.96 (0.1 x 0.9/10, 000) 1/2 = 0.1 ± 0.006. We use 
p u G {0.1, 0.5}, (7 u — 1 and values of £ suggested approximately in section 2.5 by the analyses 
of the Gulf of Mexico data (£ ~ 0.1) and the North Sea data (£ ~ —0.2). 

The plots in Figures 4 (£ = 0.1, p u = 0.5) and 5 (£ = —0.2 ,p u = 0.1) are based on simulated 
datasets of length n = 500 and n y = 10, i.e. 50 years of data with a mean of 10 observations 
per year, for N = 100,1000,10000 and 100000 years. Note that the plots on the bottom right 
have much wider y -axis scales than the other plots. It is evident that the estimative approach 
based on the MLE produces too few values in deciles 2 to 9 and too many in deciles 1 and 
10. When the true BGP distribution of Mn (from which z new is simulated in step 2) is wider 
than that inferred from data (in step 1 and using (10)) we expect a surplus of values in the 
Erst and last deciles. The estimative approach fails to take account of parameter uncertainty, 
producing distributions that tend to be too concentrated and resulting in underprediction of 
large values of z new and overprediction of small values of z new . 

The predictive approaches perform much better. Although departures from desired perfor¬ 
mance are relatively small, and vary with N in some cases, some general patterns appear. In 
Figure 4 the flat prior tends to overpredict large values and small values. The MDI prior tends 
to result in underprediction of large values. The Jeffreys prior underpredicts large values, to 
a greater extent than the MDI prior, and also tends to underpredict small values. All these 
tendencies are slightly more pronounced for £ = 0.1,p u = 0.1 (not shown). 

Figure 5 gives similar findings, although the N = 100 case behaves a little differently to the 
larger values of N. The Jeffreys’ prior is replaced by a control plot based on values sampled 
from a U(0,1) distribution. For £ = —0.2 and with small numbers of threshold excesses the 
Jeffreys’ prior occasionally produces a posterior that is also unbounded as £ £ —1/2, making 
sampling from the posterior difficult. 


10 




decile decile 

Figure 4: Proportions of simulated values of P(M N ^ z new \ x ) falling in U(0,1) deciles for 
the case £ = 0.1 and p u = 0.5. The prior is labelled on the plots. Separate lines are drawn 
for N = 100,1000,10000 and 100000. 95% tolerance limits are superimposed. 



decile 


decile 


Figure 5: Proportions of simulated values of P(M N ^ z new j x ) falling in U(0,1) deciles for 
the case £ = —0.2 and p u = 0.1. The Jeffreys’ prior is replaced by a control plot based on 
random U(0,1) samples. Separate lines are drawn for N = 100,1000,10000 and 100000. 95% 
tolerance limits are superimposed. 
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These results suggest that, in terms of predicting Mn for large N, the MDI prior performs 
better than the flat prior and the Jeffreys’ prior. However, a prior for £ that is in some sense 
intermediate between the flat prior and the MDI prior could possess better properties. To 
explore this we consider the prior (14) for 0 < a ^ 1. Letting a —> 0 produces flat prior 
for £ on the interval [—l,oo). In order to explore quickly a range of values for a we reuse 
the posterior samples based on the priors 7 Tf(ct, £) and 7 Tm(<T £)• We use the importance 
sampling ratio estimator (6) to estimate P(Mn ^ Z new \ x ) twice, once using 7 tf( 9 \ x ) as the 
importance sampling density h[6) and once using 7 tm( 6 \ x). We calculate an overall estimate 
of P(M n ^ Z new | x) using a weighted mean of the two estimates, with weights equal to the 
reciprocal of the estimated variances of the estimators (Davison, 2003, page 603). 

Figure 6 shows plots based on the MDI(0.6) prior. This value of a has been selected based 
on plots for a G (0.1, 0.2,... ,0.9}. We make no claim that this is optimal, just that it is a 
reasonable compromise between the flat and MDI priors, providing relatively good predictive 
properties for the cases we have considered. 



decile 


decile 


Figure 6: Proportions of simulated values of P(M N ^ z new j x) falling in U(0,1) deciles for 
different combinations of £ and p u under the MDI(0.6) prior. Separate lines are drawn for 
N = 100,1000,10000 and 100000. 95% tolerance limits are superimposed. 


2.5 Significant wave height data: single thresholds 

We analyse the North Sea and Gulf of Mexico storm peak significant wave heights using the 
MDI(0.6) prior suggested by the simulation study in section 2.4. We use the methodology 
proposed in section 2.1 to quantify the performance of different training thresholds. We use 
training thresholds set at the (0, 5, ..., Uk)% sample quantiles, for different Uk■ We define 
the estimated threshold weight associated with training threshold iq, assessed at validation 
threshold v(— w fc ), by 

k 

Wi(v ) = exp{T„(iq)}/ ^ exp {T v (uj)}, (15) 

i =i 

where T v (u) is defined in (8). The ratio w 2 (v)/w i(u), an estimate of a pseudo-Bayes factor 
(Geisser and Eddy, 1979), is a measure of the relative performance of threshold U 2 compared 
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to threshold u \. In section 3 these weights will be used to combine inferences from different 
training thresholds. 

The top row of Figure 7 shows plots of the estimated training weights against training 
threshold based for different iq. For the North Sea data training thresholds in the region of 
the sample 25-35% quantiles (for which the MLE of £ « —0.2) have relatively large threshold 
weight and there is little sensitivity to iq. For the Gulf of Mexico data training thresholds 
in the region of the 60-70% sample quantiles (for which the MLE of £ ~ 0.1) are suggested, 
and the threshold at which the largest weight is attained is more sensitive to iq. As expected 
from the histogram of the Gulf of Mexico data in Figure 1 training thresholds below the 25% 
sample quantile have low threshold weight. Note that for the Gulf of Mexico data the 90% 
and 95% thresholds have far fewer excesses (32 and 16) than the suggested 50. 





Figure 7: Analyses of significant wave height data by training threshold u. Top row: esti¬ 
mated threshold weights by the highest training threshold considered. Bottom row: N -year 
predictive return levels and medians of the predictive distribution of M^ for N = 100,1000 
and 10000. Left: North Sea. Right: Gulf of Mexico. 

The bottom row of Figure 7 shows that the IV-year predictive return levels and the medians 
of the predictive distribution of N -year maxima M^ are close for N = 100, where little or 
no extrapolation is required, but for N = 1000 and N = 10000 the former is much greater 
than the latter. For the North Sea data the results appear sensible and broadly consistent 
with estimates from elsewhere. From the 55% training threshold upwards, which includes 
thresholds that have high estimated training weights, estimates of the median of M 10 oo and 
Mioooo from the Gulf of Mexico data are implausibly large, e.g. 31.6m and 56.7m for the 75% 
threshold. The corresponding estimates of the predictive return levels are even less credible. 
The problem is that high posterior probability of large positive values of £, caused by high 
posterior uncertainty about £, translates into large predictive estimates of extreme quantiles. 
Figure 8 gives examples of the posterior samples of a u and £ underlying the plots in Figure 7. 
The marginal posterior distributions of £ are positively skewed, particularly so for the 95% 
training thresholds, mainly because for fixed a u , £ is bounded below by cr u /(x n — u). The 
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higher the threshold the larger the posterior uncertainty and the greater the skewness towards 
values of £ that correspond to a heavy-tailed distribution. For the Gulf of Mexico data at the 
95% threshold P(£ > 1/2 | x) « 0.20 and P(£ > 1 | cc) ~ 0.05. This issue is not peculiar to 
a Bayesian analysis: frequentist confidence intervals for £ and for extreme quantiles are also 
unrealistically wide. 

35% threshold 65% threshold 



2.0 2.5 3.0 3.5 1.0 1.5 2.0 2.5 3.0 



0.5 1.0 1.5 2.0 2.5 3.0 3.5 0 2 4 6 8 10 12 


Figure 8: Samples from 7r((x n ,£ | x), with density contours scaled to unity at the posterior 
mode (x). Dashed lines show the support of the posterior distribution: £ > a u /(x n — u ). 
Top: ‘best’ training threshold (for set at the 85% sample quantile). Bottom: 95% training 
threshold. Left: North Sea. Right: Gulf of Mexico. 

Physical considerations suggest that there is a finite upper limit to storm peak H s (Jonathan 
and Ewans, 2013). However, if there is positive posterior probability on £ ^ 0 then the im¬ 
plied distribution of H s is unbounded above and on extrapolation to a sufficiently long time 
horizon, TV) say, unealistically large values will be implied. This may not be a problem if 
Ni is greater than the time horizon of practical interest, that is, the information in the data 
is sufficient to allow extrapolation over this time horizon. If this is not the case then one 
should incorporate supplementary data (perhaps by pooling data over space as in Northrop 
and Jonathan (2011)) or prior information. Some practitioners assume that £ < 0 a priori , 
in order to ensure a finite upper limit, but such a strategy may sacrifice performance at time 
horizons of importance and produce unrealistically small estimates for the magnitudes of rare 
events. In section 3.3 we consider how one might specify a prior for £ that allows the data 
to dominate when it contains sufficient information to do so and avoids unrealistic inferences 
otherwise. 
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3 Accounting for uncertainty in threshold 


We use Bayesian model-averaging (Hoeting et al., 1999; Gelfand and Dey, 1994) to combine 
inferences based on different thresholds. Consider a set of k training thresholds Ui,... ,u^ 
and a particular validation threshold v. We view the k BGP models associated with these 
thresholds as competing models. There is evidence that one tends to get better predictive 
performance by interpolating smoothly between all models entertained as plausible a priori , 
than by choosing a single model (Hoeting et ah, 1999, section 7). Suppose that we specify 
prior probabilities P(tq),i = 1,... ,k for these models. In the absence of more specific prior 
information, and in common with Wadsworth and Tawn (2012), we use a discrete uniform 
prior P(ui) = 1 /k,i = 1 ,k. We suppose that the thresholds occur at quantiles that are 
equally spaced on the probability scale. We prefer this to equal spacing on the data scale 
because it seems more natural than an equal spacing on the data scale and retains its property 
of equal spacing under data transformation. 

Let 9i = (pj, cq, £j) be the BGP parameter vector under model tq, under which the prior is 
7i(0i | Ui ). By Bayes’ theorem, the posterior threshold weights are given by 


P v ( 'O i | x ) 


fv(. x | Uj'i ) P{Ui) 

| «j)P(U;)’ 


where /„(x \ U {) = f f v (x j 9^ tq)7r(dj | iq) d 6, is the prior predictive density of x based on val¬ 
idation threshold v under model tq. However, f v (x \ u.j) is difficult to estimate and is improper 
if n{6i | Ui) is improper. Following Geisser and Eddy (1979) we use Ilr=i fv( x r I x (r),Uj) = 
exp {T v (ui)} as a surrogate for f v (x \ u^) to give 



exp{f/(tq)} Pjuj) 
V‘_, <;xp{f,;(«,)} Pill/] 


(16) 


Let 9ij, j — 1,..., m be a sample from 7r(0,; | x), the posterior distribution of the GP pa¬ 
rameters based on threshold tq. We calculate a threshold-averaged estimate of the predictive 
distribution function of M N using 


k 

P v (M n ^ z \ x) = yP(M N ^ z | x,Ui)P v (ui | x), (17) 

1=1 

where, by analogy with (10), P(M N ^ z | x,Ui ) = (1/m) F( z 'i &ij) nyN - The solution 
zpm(N) of ^ 

P V (M 1 ^ z pm (N) | x) = 1 - 1/iV (18) 

provides a threshold-averaged estimate of the iV-year predictive return level, based on vali¬ 
dation threshold v. 

3.1 Simulation study 2: single and multiple thresholds 

We compare inferences from a single threshold to those from averaging over many thresholds, 
based on random samples simulated from three distributions, chosen to represent qualitatively 
different behaviours. With knowledge of the simulation model we should be able to choose a 
suitable single threshold, at least approximately. In practice this would not be the case and 
so it is interesting to see how well the strategies of choosing the ‘best’ threshold u* (section 
2), and of averaging inferences over different thresholds (section 3), compare to this choice 
and how the estimated weights P v (iq j x) in (16) vary over tq. 

The three distributions are now described. A (unit) exponential distribution has the prop¬ 
erty that a GP(1,0) model holds above any threshold. Therefore, choosing the lowest available 
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threshold is optimal. For a (standard) normal distribution the GP model does not hold for 
any finite threshold, the quality of a GP approximation improving slowly as the threshold 
increases. In the limit £ = 0, but at finite levels the effective shape parameter is negative 
(Wadsworth and Tawn, 2012) and one expects a relatively high threshold to be indicated. 
A uniform-GP hybrid has a constant density up to its 75% quantile and a GP density (here 
with £ = 0.1) for excesses of the 75% quantile. Thus, a GP distribution holds only above the 
75% threshold. 

In each case we simulate 1000 samples each of size 500, representing 50 years of data with 
an average of 10 observations per year. We set training thresholds at the 50%, 55%,..., 90% 
sample quantiles, so that there are 50 excesses of the (90%) validation threshold. For each 
sample, and for values of N between 100 and 10,000, we solve P v ( Mm ^ z \ x) — 1/2 
for z (see (10)) to give estimates of the median of Mm- We show results for three single 
thresholds: the threshold one might choose based on knowledge of the simulation model; 
the ‘best’ threshold u* (see section 2); and another (clearly sub-optimal) threshold chosen 
to facilitate further comparisons. We compare these estimates, and a threshold-averaged 
estimate based on (17) to the true median of Mm, 77 _1 ((1/2) 10JV ), where H is the distribution 
function of the underlying simulation model. 

The results for the exponential distribution are summarized in Figure 9. As expected, all 
strategies have negligible bias. The threshold-averaged estimates match closely the behaviour 
of the optimal strategy (the 50% threshold). The best single threshold results in slightly 
greater variability, offering less protection than threshold-averaging against estimates that are 
far from the truth. In the normal case (Figure 10) the expected underestimation is evident for 
large N: this is substantial for a 50% threshold but small for a 90% threshold. The CV-based 
strategies have greater bias than those based on a 90% threshold, because inferences from 
lower thresholds contribute, but have much smaller variability. Similar findings are evident in 
Figure 11 for the uniform-GP hybrid distribution: contributions from thresholds lower than 
the 75% quantile produce negative bias but threshold-averaging achieves lower variability 
than the optimal 75% threshold. 

In all these examples the CV-based strategies seem preferable to a poor choice of a single 
threshold, and, in a simple visual comparison of bias and variability, are not dominated clearly 
by a (practically unobtainable) optimal threshold. Using threshold-averaging to account for 
threshold uncertainty is conceptually attractive but, the exponential example aside, compared 
to the ‘best’ threshold strategy it’s reduction in variability is at the expense of slightly greater 
bias. A more definitive comparison would depend on problem-dependent losses associated 
with over- and under-estimation. 

Figure 12 summarizes how the posterior threshold weights vary with training threshold. For 
a few datasets the 90% training threshold receives highest weight. This occurs when inferences 
about £ using a 90% threshold differ from those using each lower threshold. This effect is 
diminishes if the number of excesses in the validation set is increased. In the exponential and 
hybrid cases the average weights behave as expected: decreasing in u in the exponential case, 
and peaking at approximately the 70% quantile (i.e. lower than the 75% quantile) in the 
uniform-GP case. In the exponential example the best available threshold (the 50% quantile) 
receives the highest weight with relatively high probability. In the hybrid example the 70% 
quantile receives the highest weight most often. The 75% quantile is the lowest threshold 
at which the GP model for threshold excesses is correct. The 70% quantile performs better 
than the 75% quantile by trading some model mis-specihcation bias for increased precision 
resulting from larger numbers of threshold excesses. In the normal case there is no clear-cut 
optimal threshold. This is reflected in the relative flatness of the graphs, with the average 
weights peaking at approximately the 70-80% quantile and the 50% threshold being the ‘best’ 
slightly more often than higher thresholds. Given the slow convergence in this case it may 
be that much higher thresholds should be explored, requiring much larger simulated sample 
sizes, such as those used by Wadsworth and Tawn (2012). 
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50% threshold-averaged 



Figure 9: Exponential example. Predictive median of Mm by N: individual datasets (grey), 
IV-specific (5,25,50,75,95)% sample quantiles (dashed) and true median (solid). Threshold 
strategies: median (top left); 90% quantile (bottom left); threshold-averaged (top right); 
‘best’ threshold (bottom right). 


90% threshold-averaged 



Figure 10: Normal example. Predictive median of Mm by N: individual datasets (grey), 
IV-specific (5,25,50,75,95)% sample quantiles (dashed) and true median (solid). Threshold 
strategies: 90% quantile (top left); median (bottom left); threshold-averaged (top right); 
‘best’ threshold (bottom right). 
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75% threshold-averaged 



Figure 11: Uniform-GP hybrid example. Predictive median of M N by N: individual 
datasets (grey), IV-specific (5,25,50,75,95)% sample quantiles (dashed) and true median 
(solid). Threshold strategies: 75% quantile (top left); 90% quantile (bottom left); threshold- 
averaged (top right); ‘best’ threshold (bottom right). 



Figure 12: Threshold weights by training threshold. Top: individual datasets (grey) with 
threshold-specific sample means (solid black) and (5,25,50,75,95)% sample quantiles (dashed). 
Bottom: relative frequency with which each threshold has the largest weight. Left: exponen¬ 
tial. Middle: normal. Right: uniform-GP hybrid. 
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3.2 Significant wave height data: threshold uncertainty 

We return to the significant wave height datasets, using the methodology of section 3 to 
average extreme value inferences obtained from different thresholds. Figure 13 shows the 
estimated threshold-specific predictive distribution functions of M 100 and M 10 oo- Also plot¬ 
ted are estimates from the weighted average (17) over thresholds, for different choices of the 
highest threshold Uk . For the North Sea data there is so little sensitivity to Uk that the 
black curves are indistinguishable. For the Gulf of Mexico data there is greater sensitivity 
to Uk , although based on the discussion in section 2.2 setting at the 95% sample quantile 
is probably inadvisable with only 315 observations. However, for both choices of Uk, averag¬ 
ing inferences over thresholds has provided some protection against the high probability of 
unrealistically large values of H s estimated under some individual thresholds. 



Figure 13: Threshold-specific (grey lines) and threshold-averaged (black lines) predictive dis¬ 
tribution functions of M 100 (top) and M 10 oo (bottom). Left: North Sea data. Right: Gulf of 
Mexico data. 


3.3 A weakly-informative prior for £ 

We have used prior distributions for model parameters that are constructed without reference 
to the particular problem in hand. This strategy is inadvisable when the data contain insuf¬ 
ficient information to dominate such priors, because inferences are then influenced strongly 
by a generically-chosen prior. When analysing the Gulf of Mexico data in section 2.5 we 
found that for the highest thresholds there was insufficient information in the data to avoid 
unrealistic extreme value extrapolations at long time horizons. If small sample sizes cannot 
be avoided and long time horizons are important then unrealistic inferences can be avoided by 
providing application-specific prior information. This prior could be elicited from an expert, 
as in Coles and Tawn (1996), or specified to reflect general experience of the quantity under 
study, such as the beta-type prior for £ on —1/2 ^ 1/2 used by Martins and Stedinger 

(2001) for river flows and rainfall totals. 


19 









Here, we illustrate the effects on the Gulf of Mexico analysis of providing weak information 
about £, with the aim of preventing unrealistic inferences rather than specifying strong prior 
information. The particular prior we use is 

'Kc( cr u, A) oc a” 1 (1 + £ 2 /kL 2 ) -1 a u > 0, £^ -1,A > 0, (19) 

so that the marginal prior for £ is a Cauchy distribution, truncated to £ ^ —1. The prior mode 
is zero, reflecting the expectation that £ somewhat close to zero. The use of a Cauchy density 
is motivated by Gelrnan (2006), who uses a half-Cauchy prior for a group-level standard 
deviation in a hierarchical model. The scale parameter A is set, guided by expert information 
from oceanographers (see below), to downweight a priori unrealistically large values of £. The 
tails of a Cauchy density have a gentle slope, so information from the data can be influential 
if it conflicts with this prior. 

Let tun be the median of M N . Oceanographers with expert knowledge of the hurricane- 
induced storms in the Gulf of Mexico have provided approximate values of mioo = 15m and 
1.5 for the ratio mioooo/ m ioo- Estimating mi directly from the data gives m i = 4.55m. With a 
view to downweighting substantially a apriori only very unrealistically large values of £ we set 
m 10000 = 3m 100 = 45m, that is, we double the ratio provided by the experts. Based on theory 
concerning the limiting behaviour of the maximum of i.i.d. random variables (Coles, 2001, 
chapter 3) we suppose that, for IV ^ 1, Mn has a generalized extreme value GEV(/itv, &n, £) 
distribution. Then /pv = /^i+ 04 [(IV/In2)^ — 1]/£ and the quantity R = (mioooo - mi)/(mioo — 
mi) is a function only of £. If £ = 0.229 then R = (fhioooo — mi)/(mioo — in i), suggesting a 
prior for which P(£ > 0.229) is small. We use A = 0.154, for which P(£ > 0.229) = 0.2 and 
P(£ > 1/2) = 0.1. Note that under the MDI(0.6) prior P(£ > 1/2) = 0.41. 

Figure 14 contains graphical summaries of posterior samples under the Cauchy and MDI(0.6) 
priors for 65% and 95% training thresholds. For the 65% threshold, comparing the marginal 
posterior densities for £ (top left) shows that the change of prior has little effect, because at 
this threshold there is sufficient information in the data to dominate these priors. Similarly, 
the plot of the joint posterior (bottom left) is very similar to the corresponding plot in the 
top right of Figure 8. For the 95% threshold the posterior for £ under the Cauchy prior is 
much less diffuse than under the MDI(0.6) prior and far less posterior density is associated 
with large values of £. Compare also the bottom right plots in Figures 8 and 14. Lack of 
information the data means that the posterior based on the Cauchy prior is only slightly less 
diffuse than the prior itself. 

The plots on the left of Figure 15 show, by comparison with the plots on the right of 
Figure 7, the effect of the change of prior on the threshold weights and the threshold-specific 
predictive extreme value inferences. The change from the MDI(0.6) prior to the Cauchy prior 
has changed little the threshold weights. The main difference is that the highest thresholds 
benefit more greatly from the increased prior information than the lower thresholds, resulting 
in a slight increase in their estimated threshold weights. For the highest thresholds the Cauchy 
prior has prevented the very unrealistic estimates that were obtained under the MDI(0.6) 
prior. This can also be seen by comparing the plots on the right of Figures 13 and 15: the grey 
curves corresponding to high thresholds have shifted to the left, that is, towards giving higher 
density to smaller values of M Nl with a similar knock-on effect on the threshold-averaged 
black curves. 


4 Discussion 

We have proposed new methodology for extreme value threshold selection based on a GP 
model for threshold excesses. It can be used either to inform the choice of a ‘best’ single 
threshold or to reduce sensitivity to a particular choice of threshold by averaging extremal 
inferences from several thresholds, weighting thresholds with better cross-validatory predictive 
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Figure 14: Posterior samples for the Gulf of Mexico data. Top: kernel estimates of the 
posterior density for £ under the Cauchy prior (solid) and MDI(0.6) prior (dotted), with the 
Cauchy prior density (dashed). Bottom: samples from n (cr u ,£ | x) under the Cauchy prior, 
with density contours scaled to unity at the mode (x). Dashed lines show the support of the 
posterior distribution. Left: 65% threshold. Right: 95% threshold. 


performance more heavily than those with poorer performance. The simulation study in 
section 3.1 shows that the estimated threshold weights behave as expected in cases where 
the GP model holds exactly above some threshold and illustrates the potential benefit of 
averaging different estimated tail behaviours to perform extreme value extrapolation. 

The methodology has been applied to significant wave height datasets from the northern 
North Sea and the Gulf of Mexico. For the latter dataset the highest thresholds result in 
physically unrealistic extrapolation to long future time horizons. Averaging inferences over 
different thresholds avoids basing inferences solely on one of these thresholds, but we also 
explored how the incorporation of basic prior information can be used to address this problem. 
Stronger prior information about GP model parameters, or indeed prior information about 
the threshold level itself, could also be used. 

In common with many existing methods we need a set of thresholds ui,...,Uk or an 
interval ( u m i n ,u max ) for threshold, where Uk (or u max ) does not exceed the highest threshold 
from which it is judged that meaningful inferences can be made. This is also true of other 
methods, such as those that assess formally stability of parameter estimates and those for 
which the form of a standard extreme value model is extended to model data below the 
threshold. In the latter case results can be sensitive to the form of the model specified below 
the threshold and to u min . 

The fact that our methodology is based on inferences from standard unmodified extreme 
value models makes it relatively amenable to generalization. In significant wave height exam¬ 
ples considered in this paper it is standard to extract event maxima from raw data, thereby 
producing observations that are treated as approximately independent. Otherwise, data may 
exhibit short-term temporal dependence at extreme levels, leading to clusters of extremes. 
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Figure 15: Extreme value inferences for the Gulf of Mexico data using a Cauchy prior for 
£. Top left: estimated threshold weights by the highest training threshold. Bottom left: 
IV-year predictive return levels and medians of the predictive distribution of M N for N = 
100,1000 and 10000. Right: threshold-specific (grey lines) and threshold-averaged (black 
lines) predictive distribution functions of M wo (top) and Miooo (bottom). 

In on-going work we are extending our general approach to this situation and to deal with 
other important issues: the presence of covariate effects; the choice of measurement scale; and 
inference for multivariate extremes. 
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